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Abstract 

Optimization problems associated with the interaction of linked par- 
ticles are at the heart of polymer science, protein folding and other im- 
portant problems in the physical sciences. In this review we explain how 
to recast these problems as constraint satisfaction problems such as linear 
programming, maximum satisfiability, and pseudo-boolean optimization. 
By encoding problems this way, one can leverage substantial insight and 
powerful solvers from the computer science community which studies con- 
straint programming for diverse applications such as logistics, scheduling, 
artificial intelligence, and circuit design. We demonstrate how to con- 
strain and embed lattice heteropolymer problems using several strategies. 
Each strikes a unique balance between number of constraints, complexity 
of constraints, and number of variables. In addition, each strategy has 
distinct advantages and disadvantages depending on problem size and 
available resources. Finally, we show how to reduce the locality of cou- 
plings in these energy functions so they can be realized as Hamiltonians 
on existing adiabatic quantum annealing machines. 
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1 Introduction 



1.1 Motivation and Background 

Optimization problems associated with the interaction of linked particles are 
ubiquitous in the physical sciences. For example, insights into a problem of bio- 
logical relevance such as the protein folding problem can be obtained from trying 
to solve the optimization problem of finding the lowest energy configuration of 
a given sequence of amino acids in space (Sali, Shakhnovich & Karplus 1994, 
Pande 2010, Dill, Ozkan, Shell & Weikl 2008, Mirny & Shakhnovich 2001, Pande, 
Grosberg & Tanaka 2000, Dill 1995, Gruebele & Wolynes 1998, Shakhnovich 
1994). Among other examples of biologically relevant polymers, DNA and RNA 
chains also fold into complicated structures which can be challenging to predict. 

The number of possible configurations (in fact, the number of local minima) 
for a protein with N amino acids is exponential in N (Hart & Istrail 1997). Even 
the simplest model for lattice folding (Lau & Dill 1989) was proved to be an NP- 
hard problem (Berger & Lcighton 1998, Cresccnzi, Goldman, Papadimitriou, 
Piccolboni & Yannakakis 1998). This implies that the scaling of the worst 
case scenario for arbitrary protein sequences is exponential with the size of the 
system. This scaling imposes limitations on the exhaustive search in lattice 
models for proteins with as few as 36 amino acids in even the most coarse 
grained protein models (Schram, Barkema & Bisseling 2011). 

An alternative route to exhaustive search or the development of new heuris- 
tics is to map these problems into the form of other, more general problems 
which have been extensively studied for decades. For instance, the NP-Complete 
problem known as Max-SAT has central importance to practical technologies 
such as artificial intelligence, circuit design, automated theorem proving, cryp- 
tography and electronic verification (Hansen & Jaumard 1990, Hansen, Jau- 
mard & De Aragao 1998, Soos, Nohl & Castelluccia 2009). The study of this 
particular problem is central to computer science. There are several journals, 
conferences and competitions every year dedicated entirely to solving SAT prob- 
lems (Marques-Silva & Sakallah 2007). Another widely studied constraint sat- 
isfaction problem is linear programming which has many applications including 
logistics scheduling, operations research, company management, and economic 
planning (Hemmecke, Koppe, Lee & Weismantel 2009). Some applications of 
linear programming, i.e. multi-commodity flow problems, are considered impor- 
tant enough that entire fields of research exist to develop specialized algorithms 
for their solution (Even, Itai & Shamir 1976). 

Once cast as one of these canonical constraint satisfaction problems one 
can leverage decades of progress in these fields to solve lattice heteropolymer 
problems. Though it has received relatively little attention until recently, the 
idea that constraint programming can help solve problems of this type has at 
least appeared in protein folding and computer science literature since (Yue & 
Dill 1995). Other relevant papers include (Ullah & Steinhofel 2010, Dal Palu, 
Dovier & Fogolari 2004, Krippahl & Barahona 1999, Backofen 1998, Backofen 
& Will 2006). 
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Another intriguing option is to study these problems using a computer which 
takes advantage of quantum mechanical effects to drastically reduce the time 
required to solve certain problems. For combinatorial optimization problems, 
perhaps the most intuitive quantum computing paradigms is quantum annealing 
(Finnila, Gomez, Sebenik, Stenson & Doll 1994, Kadowaki & Nishimori 1998, 
Santoro & Tosatti 2006, Das & Chakrabarti 2008, Apolloni, Cesa-Bianchi & De 
Falco 1988, Apolloni, Carvalho & De Falco 1989, Smelyanskiy, Ricffel, Knysh, 
Williams, Johnson, Thorn, Macready & Pudcnz 2012), also known as adiabatic 
quantum computation (Farhi, Goldstone, Gutmann & Sipser 2000, Farhi, Gold- 
stone, Gutmann, Lapan, Lundgren & Preda 2001, Kadowaki & Nishimori 1998). 
In quantum annealing, the presence of quantum fluctuations (tunneling) allows 
the system to efficiently traverse potential energy barriers which have a tendency 
to trap classical optimizations algorithms. 

Motivated by the experimental realization of studying biologically interest- 
ing optimization problems with quantum computation, in this contribution we 
present a general construction of the free-energy function for the two-dimensional 
lattice heteropolymer model widely used to study the dynamics of proteins. 
While the authors have already demonstrated some of these techniques in (Perdomo, 
Truncik, Tubcrt-Brohman, Rose & Aspuru-Guzik 2008), the encoding strategies 
discussed here are more general and also more efficient than what we have ex- 
plained previously. The reduction in resources achieved with these methods 
allowed for the first experimental implementation of lattice folding on a quan- 
tum device (Perdomo-Ortiz, Dickson, Drew-Brook, Rose & Aspuru-Guzik 2012) 
where we employed up to 81 superconducting qubits to solve a 6 amino- acid 
problem in the Miyazawa-Jernigan (MJ) model (Miyazawa & Jcrnigan 1996). 

The goal of this review is to explain the mapping used in (Perdomo-Ortiz 
et al. 2012), to discuss the strengths and weaknesses of this mapping with respect 
to other strategies, and to demonstrate how to map the lattice heteropolymer 
problem into forms which can be solved by using different types of technology 
and algorithms. While the focus of this paper will be on lattice protein folding, 
the methods introduced here have very general relevance to discrete and combi- 
natorial optimization problems in science. Whether one decides to use a classical 
or a quantum (annealing) device, the mappings and techniques presented here 
emphasize the importance of three key considerations: energy function locality, 
coupler/coefficient resolution, and efficiency of encoding. 

In this context, the "locality" of an expression refers to the order of the 
largest many-body expansion term. For instance, QUBO problems, which are 
a binary version of the Ising model, are said to be "2-local" because QUBO 
expressions never contain terms with more than two variables. This is a relevant 
consideration because an expression which is 3-local cannot be programmed into 
a quantum device with only pairwise couplings. A similar consideration applies 
to classical solvers. Coefficient resolution refers to the ability of a quantum 
device or classical solver to program coupler values to the degree of precision 
required for the problem. Finally, the efficiency of the encoding refers to the 
number of bits required to encode the problem. A sketch of how one might 
weigh these considerations to determine an encoding is shown in Fig. 1.1. 
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What is your computing paradigm? 



Quantum 
Annealing 



Is number of bits 
or coupler resolu- 
tion more limiting? 



Bits 



Does your device 
have programmable 
many-body terms? 



Traditional 
Computation 



Is the problem 
small enough to 
solve exactly? 




Perhaps 




Figure 1: Flow chart describing how one might choose between the three prob- 
lem encodings discussed in this review based based on available computing re- 
sources. The "diamond encoding" is not very efficient but produces a sparse 
QUBO matrix without requiring reductions that increase the required coupler 
resolution. This makes it a natural choice for classical integer-linear program- 
ming (ILP) and heuristic satisfiability (SAT) solvers which perform best on 
underconstrained problems. The "turn circuit" representation is an overcon- 
strained, but highly efficient, mapping that works best for methods designed to 
solve high-local expressions such as many-body ion trap simulators or pseudo- 
boolean optimization (PBO) solvers. The "turn ancilla" encoding represents a 
balance of these benifits as it is relatively efficient and can easily collapse to 
2-local without extremely high term coefficients. 
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1.2 Overview of Mapping Procedure 

The embedding strategies presented here apply to many discrete optimization 
problems. Mapping these problems to a constraint programming problem is a 
three step process. In this section we provide a brief description of the process 
and expand upon each step as it applies to lattice folding in later sections. 

1. Encode solution space in computational basis 

Define a one-to-one mapping between possible valid assignments of the 
problem and a bit string encoding this information. Let us denote the 
bit string by q = (?i(?2 ,,, 9n- The way information is encoded at this 
point can drastically alter the nature of the following three steps so one 
must take care to choose a mapping which will ultimately make the best 
use of resources; in many cases, the most compact mapping will have a 
high order energy function or require many ancillary bits. Regardless of 
how information is encoded, the bit string must uniquely enumerate each 
element of the low energy solution space. 

2. Constrain energy landscape with pseudo-boolean expression 

Construct a pseudo-boolean energy function E(q) = E(q\,q2,--- ,q n ) 
which takes q as input and correctly reproduces the relative energies in the 
low energy subspace of the original problem so that the optimal solution 
to E{q) encodes the solution to the original problem. The construction 
of this function is not trivial and will depend largely on how informa- 
tion is encoded in q. At this point it may be necessary to increase the 
dimensionality of the solution space by adding ancillary bits. In a previ- 
ous contribution, we provided a specific technique to construct the energy 
function for particles interacting in a lattice (Pcrdomo et al. 2008). The 
purpose of this contribution is to introduce the reader to several differ- 
ent types of mappings which have distinct advantages or disadvantages 
depending on problem size, complexity and available resources. 

3. Map boolean representation to desired constraint programming 

In most cases one can take advantage of significantly more powerful solvers 
by making a final transformation from pseudo-boolean function to weighted 
maximum satisfiability (W-SAT), integer-linear programming (LLP), or 
quadratic unconstrained binary optimization (QUBO). When cast as a W- 
SAT problem one can take advantage of both heuristic and exact W-SAT 
solvers which have been developed by the computer science community 
and tested every year in annual "SAT Competitions" . When represented 
as an ILP problem, one can use commercial logistics scheduling software 
such as IBM's CPLEX. If one wishes to implement the energy expression 
on a quantum device it may be necessary to manipulate the energy ex- 
pression so that it contains only local fields and two-body couplings. So 
the final step is often to reduce the dimensionality of the pseudo-boolean 
expression to 2-local so that the problem can be implemented as QUBO 
on currently existing architectures for adiabatic quantum computing as 
was done in (Perdomo-Ortiz et al. 2012). 
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2 The "Turn" Encoding of Self- Avoiding Walks 



2.1 Embedding physical structure 

Let us use the term "fold" to denote a particular self-avoiding walk (SAW) 
assumed by the ordered chain of beads or "amino acids" on a square lattice. 
These configurations include amino acid chains that might intersect at differ- 
ent points due to amino acids occupying the same lattice sites. Even though 
overlapping folds will exist in the solution space of our problem, these folds are 
unphysical and therefore we need to construct energy functions to penalize such 
configurations. Such functions will be discussed in detail below. 

A fold of an N amino acid protein is represented in what we refer to as the 
"turn" mapping by a series of A — 1 turns. We use this name to distinguish the 
encoding from other (spatial) representations which encode the possible folds 
by explicitly encoding the grid location of each amino acid. The square lattice 
spatial representation discussed in (Perdomo et al. 2008) has the advantage of 
being general for the problem of N particles interacting in a lattice (which need 
not be connected) but we can do much better in terms of the number of variables 
needed; bit efficiency is the main advantage of the turn mapping. 

In the turn mapping, one saves bits by taking advantage of the connectivity 
of a valid SAW to store information about where each amino acid is relative 
to the previous amino acid instead of encoding explicit amino acid locations. 
Therefore, instead of encoding the positions of the jth amino acids in the lattice, 
we encode the jth turn taken by the j + 1 amino acid in the chain. For peda- 
gogical purposes, we concentrate on the case of a two-dimensional (2D) lattice 
SAW; the extension to a three-dimensional lattice requires a straightforward 
extension of the same techniques described here for the 2D case. 

Because the location of an amino acid in the turn mapping is specified by 
its location relative to the previous acid in the primary sequence, the solution 
space consists only of paths, or "worms" , embedded in the lattice. The result- 
ing energy function is translationally and rotationally invariant with respect to 
the embedding in physical space as long as the local structure of the relative 
locations is kept intact. More specifically, each of the N—l turns in 2D space re- 
quires two bits so that each of the four directions (up, down, left, and right) has 
a unique representation. This assumes a rectilinear lattice, but the method is 
equally valid, though with slight modification, for other lattices, e.g. triangular. 
The convention or "compass" used in this paper is presented in the upper-left 
part of Fig. 2. Furthermore, we can fix the first three bits to obtain only solu- 
tions which are rotationally invariant. Under this convention, the bit-string q 
is written as, 



We have chosen to fix the first three bits as 010 so that the walk always 
turns first to the right and then either right or down. This does not affect 
the structure of the solution space and leaves only N — 2 turns to be specified; 




(1) 



turn2 tumZ 



turn(N-l) 
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an example is provided in Eq. 1. Since every turn requires 2 bits, the turn 
mapping requires only 2(N — 2) — 1 = 2N — 5 bits to represent a fold. This 
can be compared with the (2N — 4) log 2 N required for the spatial mapping in 
(Perdomo et al. 2008). To clearly demonstrate how this mapping works, an 
example of the turn encoding for a short SAW is shown below in Fig. 2. 
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Figure 2: Step-by-step construction of the binary representation of a particular 
six unit lattice protein in the turn encoding. Two qubits per bond are needed 
and the turn "compass" (bond directions) are denoted as "00" (downwards), 
"01" (rightwards), "10" (left), and "11" (upwards). This image has been repro- 
duced from (Perdomo-Ortiz et al. 2012) with permission from the authors. 



2.2 "Turn ancilla" construction of E(q) 

Now that we have a mapping function which translates a length 2N — 5 bit- 
string into a specific fold in the 2D lattice we can construct E(q) as a function 
of these binary variables. For the case of lattice folding, we need to penalize 
folds where two amino acids overlap, i.e. the chain must be self-avoiding. This 
penalty will be given by the energy function, E over i ap {q), which returns an 
extremely high value if and only if amino acids overlap. While it is possible to 
construct a single function E overlap (q) which penalizes all potential overlaps, 
we will show that less ancillary bits are needed if we introduce the function 
Eback (q) which penalizes the special case of overlaps that happen because the 
chain went directly backwards on itself. In this scheme, E over i ap (q) will apply 
to all other potential overlaps. 
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Finally, we must consider the interaction energy among the different amino 
acids. This will ultimately determine the structure of the lowest energy fold. 
The energy given by the pairwise interaction of beads in our chain will be given 
by E pair (q). In some lattice protein models such as the Hydrophobic-Polar 
(HP) protein folding model, there is only one stabilizing interaction; however, 
the construction we present here applies for an arbitrary interaction matrix 
among the different amino acids (or particles to be even more general). One of 
the advantages of the turn representation over the spatial representation is that 
we do not need to worry about having the amino acids linked in the right order 
(primary sequence), since this is guaranteed by design. The construction of the 
energy function, 

E(q) = E back (q) + E overlap (q) + E pair (q), (2) 
involves a series of intermediate steps which we outline next. 

2.2.1 Construction of Eb ac k(q) 

In order to have a valid SAW we need to guarantee that our "worm" does not 
turn left and then immediately turn right or vice versa or turn up and then 
immediately turn down or vice versa. In order to program this constraint into 
the energy function we will introduce several simple logic circuits. Looking at 
the compass provided in Fig. 2 it should be clear the circuits in Figs. 3-6 return 
true if and only if a particular turn (encoded (71(72) went right, left, up, or down 
respectively. 



Figure 3: A logical circuit representing "right" consisting of a not gate after 
the first bit and an and gate. Evaluates to true if and only if qi, q 2 = 0, 1. 

9l 92 



(92 -9192) 




Figure 4: A logical circuit representing "left" consisting of a not gate after the 
second bit and an AND gate. Evaluates to true if and only if (71, q 2 — 1,0. 
9i 92 



(91 - 9192) 
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Figure 5: A logical circuit representing "up". Only true if q\, q% = 1, 1. 
9i 92 




(?i92) 



Figure 6: A logical circuit representing "down". Only true ii qi,q 2 = 0,0. 

qi 92 



V 




(1 - 91 - 92 + 9192) 



Using these circuits we can generalize the concept of "up", "down", "left" 
and "right" functions to precise directional strings. In two dimensions (as pre- 
scribed by Fig. 2), we have the functions for the jth turn, 

d J x+ = 92j(l-92j-i) = 92j-92j92j-i (3) 

4- = (1 - = Q23-1 ~ 92j92j-i (4) 

4+ = 92j92i-i (5) 

d J v - = (l-92j)(l-92j-i) = l-92j-92j-i+92j92j-i, (6) 

which evaluate to true if and only if the jth turn is to be right, left, up 
or down respectively. Having defined these circuits we can construct a more 
complicated circuit which takes two turns (the 4 bits 9i9i+i9i+29i+3) as in- 
put and returns true if and only if the second turn went backwards, i.e. 
(d{ + A V (d{_ A V (d? y+ A rf^t 1 ) V (d? y _ A e^ 1 ) . An example of 

these conjunctions, (d J x+ A d^X 1 ^ is shown in Fig. 7. 

The other three conjunctions are also trivially constructed by combining 
the appropriate circuits using and gates which simply multiply together the 
directional strings. The utility of these circuits is that they produce terms in a 
pseudo-boolean function. Specifically we get the terms, 

( d i+ A 4- 1 ) = - g l «?i+igi+2 - gi+ig l +29 l +3 + <?i? l +i9i+2«?i+3 (7) 

(d 3 x _ A diX 1 ') = ^+3-^+1^+3-^+2^+3 + ^+1^+2^+3 (8) 

( d a+ Ad y- 1 ) = 9»9i+i - 9i9i+i9<+2 - 9»9i+i9i+3 + 9i9i+i9<+2gi+3 (9) 

( dJ y- Ad l+ 1 ) = «i+2?i+3 - 9i?i+2?i+3 - 9i+l9i+2?i+3 + ?i9«+ig»+2gi+3- (10) 

It might seem logical to finish this circuit by combining all four backwards 
overlap circuits with OR gates; however, this is not an advisable strategy as it 
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<li Qi+2 <7i+3 



V 






Figure 7: A logical circuit which returns true if and only if (d J x + A d^t 1 ^, i.e. 
the turn sequence Qiqi+iqi+2qi+3 = 0110, meaning it went right and then left. 



is sure to produce many high ordered terms. Because exactly one or none of 
these circuits will be true we can accomplish the same result by summing the 
four circuits. Accordingly, for the two turns <Mi+i<7i+2?i+3 the pseudo-boolean 
expression, 

(di + A dit 1 ) + (4- A 4?) + (<+ A djl 1 ) + (<_ A (11) 

evaluates to true if and only if qtqi+iqi+2qi+3 represents a backwards turn and 
evaluates to false otherwise. Our goal is to construct a pseudo-boolean ex- 
pression which returns a penalty whenever a backwards turn is made; therefore 
we must multiply this expression by a constant to be determined later known 
as X over i a p. After substituting Eqs. 7-10 into Eq. 11, factoring the terms, and 
adding in \ ove ria P we can write, 

Eback {qiqi + iqi + 2qi+3) = ^overlap { < 2-qiq%+2 — qi — qi + 2) (2^ + 1^+3 — qi+i — qt + s) ■ 

(12) 

To construct the entire -E& ac fc (q) we need to sum together bits from each 
pair of adjacent turns. Keeping in mind that we fix the first three bits at 010, 
we write the final expression for Eb ac k{q) as, 



E back (q) = Xback (qi qi + qiqz - 1q\qiqz) 

+ ^back Si=2 (2giQi+2 — qi — qi+2) (2qi+iqi+3 — qi+1 — qi+s) ■ 



(13) 



In this expression, the first three terms come from ensuring that the second turn 
(which begins with a bit fixed at 0) does not overlap with the third turn. Notice 
that in this expression, the first physical bit with an unknown value is labeled 
u qi" despite the fact that the first three information bits are fixed at 010. This 
formalism will be consistent throughout our review. 

It is important to point out that while the decision to use a separate Eb ack (q) 
instead of a more general E over i ap (q) has the disadvantage of introducing 3 and 
4-local terms, it has the advantage of construction without any ancillary bits. 
Furthermore, even if one needs an entirely 2-local expression this strategy may 
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still be preferable because the same reductions needed to collapse this expression 
to 2-local will be needed in collapsing the pairwise energy function to 2-local by 
construction. For more on reductions, see Sec. 6. 

2.2.2 Construction of E over i ap (q) with ancilla variables 

The overlap energy function E over i ap (q) penalizes configurations in which any 
two amino acids share the same lattice point. The penalty energy associated 
with any pair of amino acids overlapping must be large enough to guarantee 
that it does not interfere with the spectrum of the valid configurations (we 
return to the topic of choosing penalty values later on). We begin by defining a 
function which specifies the x and y grid positions of each amino acid. Because 
the directional strings we defined earlier in Eqs. 7-10 keep track of the direction 
of every step we can define these functions as, 

n-1 

X n = l + qi + J2(4+-€-) (14) 
k=2 
n-1 

y n = qi -l + J2(d k v+ -d k y _) (15) 

k=2 

where the position of the nth amino acid in the sequence is a function of the 
proceeding n — 1 turns iterated through with index k. Note that the terms 
in front of the sum are determined by the first three (fixed) bits: 010. With 
these definitions we can make an extremely useful function which will return 
the square of the grid distance between any two amino acids (denoted i and j): 

9ij = {xi- Xjf + (yi- yjf . (16) 

gij has several extremely useful properties worth pointing out now. First, gij 
is zero if and only if two amino acids overlap; otherwise, g^ is always positive. 
Additionally, has the very surprising property of being natively 2-local when 
constructed using the compass that we defined in Fig. ?? (therefore the deci- 
sion to encode directions in that fashion was not arbitrary). This is surprising 
because the directional strings are 2-local so we might naively expect something 
which involves the square of these to be 4-local; however this turns out not to 
be the case because x n and y n are 1-local by construction. 

In order to use g^ to construct E over i ap (q) we need a function which takes 
gij as input and returns a penalty if and only if g^ = 0. First, we note the 
bounds on g^ , 

0- • (t jr- (17) 

To help enforce the constraint that g^ > 1, we introduce a free parameter, a^ . 
In the optimization literature, such a variable is called a "slack variable" and is 
used to convert an inequality into an equality. In our case, 

0<ay<(i-j) 2 -l (18) 
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This implies that, 

Vftj > 1 3 atij : (i - j) 2 - .gy - a tj = 0. (19) 
Furthermore, if and only if gy = 0, 

(*-j) 2 -ftj-a«>lVa . (20) 

In order to introduce a slack variable such as ay into the construction of our 
pseudo-boolean function we must encode it using ancilla bits. Ancilla bits are 
real, unconstrained bits used in the calculation which have no physical signifi- 
cance to the particular problem mapping (i.e. ancilla bits do not tell us anything 
about a particular protein fold) . In using ancilla we increase the dimensionality 
of the solution space of our problem by introducing extra variables but gain the 
ability to use those bits in our energy function. 

Every pair of amino acids which could possibly overlap will need unique 
bits to form an a for use in the E over i ap (q) term corresponding to that pair. 
Only amino acids which are an even number of turns apart can possibly overlap 
and we are already preventing amino acids which are two turns apart from 
overlapping with Eb ac k(q); thus, the number of amino acid pairs which require 
a slack variable is calculated as, 

N-4 N 

E E [(1 mod 2]. (21) 

i=l j=i+4 

Each ay can be represented in binary using the corresponding ancilla bits. 
Using Eq. 18 we see that the ay corresponding to amino acid pair i,j can be 
represented in /zy ancilla bits where, 

Hij = \2 log 2 (i - j)] [(1 + i - j) mod 2] . (22) 

Therefore, the total number of ancilla bits required to form E over i ap (q) is, 

N-4 N 

E E wj- ( 23 ) 

i=l j=i+4 

Finally, we can write the formula for a given ay as, 

a« = E«c«+fc2 fc (24) 

fe=0 

where cy denotes a pointer to the first ancilla bit corresponding to a particular 
amino acid pair. For instance, if the E over i ap (q) ancilla are in sequential order 
from lowest index pair to highest index pair and come immediately after the 
information, bits then we could write, 

i / N \ N 

c ij = ^ ' I ^ ' Cmn J ~ ^ ' Min- (25) 
m— 1 \n— m+4 / n—j 
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However, there are still several problems we must address before we can 
construct E over i ap (q). To begin with, we originally wanted an djj which was 
specifically restricted to the domain given in Eq. 18 but since we cannot con- 
strain the physical bits in any fashion, Eq. 22 and Eq. 24 suggest that our slack 
variable is actually in the domain given by, 

< ay < 2^ - 1. (26) 

We should adjust Eq. 19 and Eq. 20 so that, 

Vgij > 1 3 ay : 2^ - 9ij ~ ay = 0. (27) 

Furthermore, if and only if <7y = 0, 

2"«-ffy-a tf >lVa tf . (28) 

Finally, there is the question of how to guarantee that ay is the particular ay 
that gives in Eq. 27 whenever gy > 1. Even though there exist ay such that 
Eq. 27 evaluates to 0, it is also possible to have ay such that Eq. 27 evaluates to a 
negative value. Negative values would incentivize overlaps instead of penalizing 
them so to ensure that the lowest energy solution always has E over i ap (q) = 
we square the expression to obtain the following formula, 

7y — ^overlap [2^ J Qij ^ij] • (^^) 

The expression 7y is effective for our purposes because ay's restricted domain 
given by Eq. 26, promises that 7y can only equal zero if #y > 1. 7y is zero 
only if gy > 1 A ay = 2^ ij — gij\ thus, the goal is to make Xoveriap a sufficiently 
large penalty that all low energy solutions must have no overlaps, i.e. <7y > 1 
for all ij, and ay = 2^ ij — g^. Finally we can write the final expression, 

JV-4 N 

E ove rla P ^ + * ~ ^ m ° d ^ ^ ' ( 30 ) 

i=l j=i+i 

Again, we include the term [(1 + i — j) mod 2] because only amino acids that 
are an even number apart have the possibility of overlapping. Furthermore, 
because overlaps between adjacent amino acids are impossible and overlaps be- 
tween amino acids two apart are prevented by Ef, a ck (q), we start the second 
sum at j = i + 4 Accordingly, one should only create ancillary bits for pairs in 
which (i — j) mod 2 = A i — j > 4. It should now be clear that the reason we 
introduced Eb ac k (q) was so that we used fewer ancillary bits in this step. 

2.2.3 Construction of E pair (q) with ancilla variables 

Finally, we need to construct the pairwise interaction energy function. To do 
this we need to make an interaction matrix, J, which contains all of the pair- 
wise interactions which lower the energy when two amino acids are adjacent 
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on the lattice (thus all elements of J are negative or zero). Note that this in- 
teraction matrix must contain many zero-valued elements as many amino acid 
pairs cannot possibly be adjacent. For instance, only amino acids which are at 
least three turns apart and an odd number of turns apart can ever be adjacent. 
Furthermore, depending on the interaction model many of these amino acids 
might not "interact" ; for instance, in the HP-model only H-H pairs can interact 
where as in the Miyazawa-Jernigan model all amino acids can interact. 

For each potential interaction, we must introduce one ancillary bit denoted 
Uij where i and j denote the amino acids involved in the interaction. Wjj is 
essentially a switch which is only "on" without incurring an energy penalty if 
two amino acids are interacting (that is, if g^ — 1). We can now write the 
pairwise interaction term: 

Vij = ^ijJij (2 - g^) (31) 

This simple function does everything we need to write the pair function. 
Because E overlap (q) ensures that > 1, we see that (fij is only positive if 
both Jij and uiij are non-zero and gij is greater than 2. Such solutions will 
never be part of the low-energy landscape for our problem because the energy 
could be made lower by trivially flipping the uJij ancillary bit. On the other 
hand, ifij = Jij if and only if g^ = 1 A u>ij = 1 which means that the pair is 
adjacent! Thus, the final form of E pair (q) is, 

N-l N 

Epair(q) = ^2 ^2 J i:j (2 - g^) . (32) 

i=\ j=i+3 

2.3 "Turn circuit" construction of E(q) 

The turn ancilla construction has the advantage of providing an energy expres- 
sion with relatively few many-body terms but it does so at the cost of intro- 
ducing ancilla bits. If one intends to use a pseudo-boolean solver or a quantum 
device with adjustable many-body couplings, bit efficiency is much more im- 
portant than the particular structure of the energy expression. This section 
explains the so-called "circuit" construction which provides optimal efficiency 
at the cost of introducing high ordered many-body terms. The turn circuit 
construction (along with reductions explained in Sec. 6) was used to encode 
problems into a quantum annealing machine in (Perdomo-Ortiz et al. 2012). 

2.3.1 Sum strings 

The circuit construction works by keeping track of the turns in between amino 
acids to determine if the amino acids overlap or not. To do this we keep track 
of the turns in every direction using the directional strings defined in Eqs. 7-10. 
Using these directional strings we introduce ancillary bits referred to as "sum 
strings". Sum strings are strings of [log 2 (j — i)~\ bits for each segment of the 
chain between amino acids i and j, with 1 < i < j < N and i + 1 < j. As in 
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the case of the directional strings, we require one 'sum string" per direction per 
pair of amino acids to be compared. Each represents, in binary, the number of 
total turns in a particular direction within the segment. 

As in the ancilla construction, whether or not two amino acids interact or 
overlap depends on the sequence of turns between them. To determine this, for 
each segment of the directional strings we construct a string that is the sum, in 
binary, of the bits between two amino acids, i.e. the total number of turns in 
that direction. This process is most straightforwardly described using a circuit 
model. Consider, a single Half- Adder gate (HA) consisting of an and and an 
XOR gate, as shown in Fig. 8. The output of a Half-Adder can be interpreted 



x - 

y—\ 



HA 



- x Ay 

— x (By 



Figure 8: The Half- Adder gate sums two bits. 



as the two-bit sum of its two input bits. Accordingly, if we wanted to add three 
bits we could add two of them, and then add the resultant two-bit number to 
the third bit, as shown in Figure 9. 



HA 



- x Ay - 
-x®y- 



HA 



- (x © y) A z - 



HA 



-0 



(lAj)e {(x ®y)f\z) 

■ (x © y) © z = (x + y + z) mod 2 



Figure 9: A circuit to sum three bits. 



In general, to add a single bit to an n-bit number, we simply apply n Half- 
Adders. First, a Half-Adder applied to the single bit and the least significant 
bit of the augend gives the least significant bit of the sum. Next, we use a 
second Half- Adder to add the carry bit of the first addition and the second least 
significant bit of the augend to give the second least significant bit of the sum. 
This process is repeated until the (n + l)-bit sum is computed. For an example 
of this see Fig. 10. 
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x 5 ■ 
X4 - 

^3 - 
X2 ■ 
X\ - 

y — I 



HA 



HA 



HA 



HA 



HA 



\—z e 

■ z 5 

■ z 4 

Z3 

■ z 2 

■ Zl 



Figure 10: Circuit for the addition of a single bit y to the 5-bit x = x 5 X4X^X2Xi 
to form the 6-bit sum x + y = z 6 z 5 Z4Z 3 z 2 zi. 

Thus, given an arbitrary number of bits we can find their sum, in binary, 
by successively combining the strategies shown in Fig. 11, i.e. first adding the 
first three bits (see the first three HA gates from left to right) and then adding 
the next bit to the resulting three bit number which carries the previous sum. 
This is accomplished by the next three HA gates. From then and on, one adds 
a simple bit to each of the resulting n bit number by using n HA gates until all 
bits in the string are added. 



x 5 



II A 



£3 

Xl 



HA 



HA 



HA 



HA 



HA 



HA 
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HA 
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■ s 4 
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Figure 11: The circuit for the sum, S1S2S3S4S5, of 5 bits x\ + x 2 + x 3 + x 4 + x 5 . 
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-4,±(m) 
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Figure 12: The circuit for the number Sk±(i,j) of turns between amino acids i 
and j in the ±fc direction. 



We can use the circuit such shown in Fig. 12 to compute the the binary 
digits of a particular sum that will be very useful to us, 



s r k± (i,j) = r th digit of 



(33) 



This sum tells us how many turns our protein has taken in the k± direction 
between any two amino acids. For instance, s*_ (3, 9) would tell us the value 
of the 1st binary digit of an integer representing the number of times that the 
protein turned in the negative x direction (aka left) between amino acids 3 and 
9. While the size of the output of the circuit given in Fig. 12 scales exactly with 
the size of the input, the maximum number of bits needed to represent the sum 
of a set of bits scales logarithmically; therefore, many of the bits representing 
higher places in the sequence arc zero. Specifically, the sum of n bits requires 
at most [~log 2 n] bits to represent in binary. 



2.3.2 Construction of E over i ap (q) with circuit 

The overlap penalty should be positive if any two amino acids are at the same 
lattice point. For a pair i,j, this occurs when the number of turns between 
them in each direction fc± is equal to those in the opposite direction fc=p or 
equivalently, when the bit-strings representing those numbers, s J k ~^ ■ ■ ■ sj. + and 
sjT 1 ■ ■ ■ s l-, are the same. As discussed above, since only the first [log 2 (j — 
digits of Sk± are non-zero, the overlap penalty function for amino acids i,j is 

D /[log(j-i)l \ 

E overlap ^j) = [] II XNOR(4 + (i,j),sUU)) , (34) 

fe=l y r=l J 
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where 



XNOR(p, q) = l-p-q + 2pq 



(35) 



is the exclusive NOR function which evaluates to true if and only if the two bits 
have the same value. Furthermore, we need not consider every pair of amino 
acids in the sequence because in order for the number of turns in opposite 
directions to be equal, there must be an even number of total turns. The total 
on-site penalty function is 

N-2 /L(JV-*)/2j \ 
^overlap — ^overlap ^ ^ J ^ ^ ^overlapi^i % H~ 2j) I (36) 

»=i \ j=i / 

2.3.3 Construction of E pair (q) with circuit 

To determine if a pair of amino acids is adjacent on the lattice without using 
ancilla bits is more involved. Two amino acids are adjacent if and only if the 
number of turns between them in opposite directions is the same in all but one 
dimension and the numbers of turns in the other dimension have a difference 
of one. The construction of the equality condition is the same as in as for the 
overlap function; to construct the latter condition, consider the set of 4 bit 
numbers, as shown in Figure 13. 



Decimal 


Binary 





0000 


1 


0001 


2 


0010 


3 


0011 


4 


0100 


5 


0101 


6 


0110 


7 


0111 


8 


1000 


9 


1001 


10 


1010 


11 


1011 


12 


1100 


13 


1101 


14 


1110 


15 


1111 



Figure 13: All 4 bit binary numbers and their decimal representations. 

Note that when the first of two sequential binary is even, the Hamming 
distance between those bit-strings are the same except for the least significant 
bit, e.g. 0000 and 0001, 1000 and 1001, 1110 and 1111. On the other hand, 
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sequential numbers for which the lesser one is odd differ in at least two places, 
depending on where the rightmost is in the lesser number, i.e. 



00000000000 • • • 01 
+ **...* *011 • • • 11 , (37) 
**---**100---00 

as in 0011 and 0100, 0111 and 1000, and 1011 and 1100. 

Let us use p to denote the position of the rightmost in the odd, lesser 
number of this comparison. There are three portions of the bit strings which 
need attention when comparing adjacency in this case. First, all digits from 
the least significant and up to p need to be different. Second, all digits after p 
need to be the same. Third, within each possible adjacency direction (k+ or 
k— ) there needs to be a change from p — 1 to p. Finally, all the digits from 
the least significant up to the (p — 2)th digit need to be the same. Using these 
conditions, for both cases when the lesser number is cither even or odd, results 
in the adjacency terms for each of the two dimensions and all of the possible 
amino acid pairs, afe(i,j): 



ak(i,j) = 



I riogo-oi N 

n n xnor (s r w+ (i,j), s r w _(i,j)) 



(38) 



XOR(4+(U),4- EI XNOR(4+(;,j),4-(i,j)) 



r = 2 

ri°g(j-*)i / p-2 



+ 



3(3-!) I / P-2 

Y, XOR^; 1 ^),^^)) nXNOR^+^j),^ 1 ^)) 



p = 2 \ r=l 
P [log(j-i)] N 

*UX0R(s r k+ (i,j),sl_(i,j)) XNOR (sl + (i,j),sl_(i,j)) 

r— 1 r— p+1 y 

Thus total contribution of the interaction between two amino acids to the 
total energy function is given by 

E P air(i,j) = Jij [ax(i,j) + a y (i,j)] , (39) 

where Jij is the adjacency matrix giving the energy of pairwise interactions that 
we used earlier. As was the case with the overlap penalty function, we need not 
consider all pairs of amino acids. In order for two amino acids to be adjacent 
there must be an odd number of turns between them, excluding the trivial case 
of amino acids that are adjacent in the primary sequence. Accordingly, the total 
pairwise interaction function is 



E pair = V V E pair (t, l + i + 2j) ) . (40) 
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3 The "Diamond" Encoding of SAWs 



There are many different ways in which one couid encode a SAW into binary. Of 
all the alternatives to the "turn" encoding that we have considered, one stands 
out for a number of reasons: the so-called "diamond encoding" lends itself 
to an energy function which is natively 2-local (without any reductions) and 
which has a very sparse QUBO (quadratic unconstrained binary optimization) 
matrix. Despite the fact that the diamond encoding requires no ancillary bits 
whatsoever, the encoding is still less bit-wise efficient than the "turn encoding" . 
In the language of constraint satisfaction programming, this means that the 
clause:variable ratio is significantly lower when compared to the clause:variable 
in the turn encoding. 

3.1 Embedding physical structure 

The diamond encoding can be thought of as a more sophisticated and subtle 
version of the "spatial" encoding used in (Perdomo et al. 2008). The key insight 
behind the diamond encoding is that if the first amino acid is fixed then each 
subsequent amino can only occupy a very restricted set of lattice points which 
can be enumerated independent of any knowledge of the particular fold. To 
clarify this point and elucidate why we refer to this as the "diamond" encoding, 
see Fig. 14. 




Figure 14: A "map" of the diamond encoding in 2D. If the first amino acid is 
fixed to the blue lattice point in the center then the second amino acid must be 
on an orange lattice point, the third must be on a green lattice point and the 
fourth must be on either an orange or red lattice point. 

Fig. 14 illustrates what the "diamond" of valid lattice sites looks like for 
the first 4 amino acids in a SAW. In the diamond encoding each bit refers to 
a specific lattice site which could be occupied by an amino acid in that part of 
the sequence. In Fig. 14 we notice that the second amino acid may occupy 4 
positions, the third may occupy 8 and the fourth may occupy 18. Accordingly, 
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we need this many bits for each amino acid. 

9 = gig2g3g4 g5g6£rg8g9£lQgllgl2 ■ ■ ■ (41) 
2nd acid 3rd acid 

Though very straightforward to encode, this representation makes signifi- 
cantly less efficient use of bits than does the turn representation. However, 
there are a few tricks which we can use to improve the situation for this encod- 
ing. While the "diamond" of possible lattice locations for each amino acid grows 
quadratically with the length of the chain we can simultaneously save bits and 
drastically reduce the solution space without discarding the global minimum by 
deciding to set a hard limit on the size of the diamond. For instance, if a protein 
has length 100 then we would expect that the diamond for the 100th amino acid 
will have a radius of 99 lattice points at each corner. However, we can use the 
observation that proteins always fold into very compact structures to justify a 
substantial restriction on the solution space of our problem. 

The very fact that these problems are typically called "protein folding" sug- 
gests that low energy solutions involve dense conformations. Indeed, almost 
all heuristic methods for folding proteins take advantage of the compact na- 
ture of low energy folds to constrain search procedures (Baker 2000, Oakley, 
Wales & Johnston 2011, Shakhnovich 1996). A large part of the reason why lat- 
tice heteropolymer problems such as protein folding are so difficult and poorly 
suited to heuristic algorithms is because the low energy solutions are always 
very compact and thus, frustrated, which makes it very unlikely that compact 
folds will be found efficiently via stochastic searches (Dill 1993, Shea, Onuchic 
& Brooks 2000, Camacho 1995, Nymeyer, Garcia & Onuchic 1998). Therefore, 
for any interesting problem its reasonable to assume that the protein will not 
stretch out further than a certain limit. To estimate this limit one must have 
familiarity with the types of solutions expected of the particular problem. An 
examination of several publications holding current records for lowest energy 
folds in canonical problems suggests to us that for a 100 unit instance in 2D a 
reasonable cutoff radius would be around 20-30 lattice points. The cutoff ra- 
dius could reasonably be made shorter for lattice models in higher dimensions 
as folds are expected to be even more compact on higher dimensional lattices. 
The number of bits required for the diamond encoding can be expected to grow 
cubicly up to a limit and then linearly after that limit if a cutoff is imposed. 
Because the number of bits required for the turn ancilla grows quadratically, for 
large proteins or proteins on higher dimensional lattices the diamond encoding 
would actually be more efficient in bit resources. 

3.2 Natively 2-local E(q) 

The major advantages of the diamond encoding become evident as soon as one 
starts to construct E(q). The breakdown of the energy function looks different 
for the diamond encoding than for the turn encoding because the diamond 
encoding has different strengths and weaknesses. The first difference is that the 
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diamond encoding will require a constraint, E one (q) which makes sure that each 
amino acid will have only one bit flipped to "on" so that each amino acid can 
only occupy one lattice position. Furthermore, the diamond encoding does not 
hard-code a primary structure constraint so we will need a term, E connect (q) to 
guarantee that each sequential amino acid is adjacent. Like the turn encoding 
the diamond encoding will also require E over i ap (q) and E pair (q) terms. The 
overall energy function looks like, 

E (q) — E one (q) + E connect (q) + E over i ap (q) + E palr (q) . (42) 

3.2.1 Construction of E one (q) 

Each amino acid is encoded by flipping a bit in the part of the total bit-string 
sequence which represents that amino acid. Thus, we need to make sure that 
exactly one bit is flipped "on" for each amino acid. The most efficient way 
to guarantee this is the case for low energy solutions is to lower the energy 
whenever a bit is flipped on but introduce extremely high penalties if any two 
are flipped on for the same amino acid. For instance, if q k is the binary vector 
which represents the fcth amino acid and nu represents the length of this vector 
then we can write, 

N n k -l n k 

E one (q) = x one EEE «M ■ ( 43 ) 

k=2 »=1 j>i 

\ one in Eq. 43 yields terms which impose a very large penalty if any two (or 
more) bits are flipped at once. As written, this function allows for the possibility 
that no bits are flipped on at once (and clearly one must be flipped on). However, 
the terms introduced in E connect (q) will guarantee that the low energy solutions 
all have one bit flipped on. Thus, this function only needs to make sure that no 
more than one bit is flipped for each amino acid. 

3.2.2 Construction of E connect (q) 

To form E connect (q) we take a very similar approach to how we formed E one (q). 
To guarantee that the low energy solution space contains only amino acids chains 
which connect in the desired order we couple every bit representing amino acid 
k to each of the n^-i < 4 bits representing a lattice position adjacent to that 
amino acid from the previous amino acid k — 1 and multiply by a reward as 
follows (using the same notation as was used in Eq. 43, 

N Tifc-lnfc-i 

Econnect (q) = N — 2 — \ CO nneet E E E ^ ^ ' ( 44 ) 

k=2 i=l j=l 

Note a subtle difference between the second and third sums here is that the 

in the upper limit of the sum is subscripted in the latter but not in the former 

equation. Another important caveat is that A co „„ ect << A one so that the system 
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cannot overcome the \ one penalty by having multiple A cfmnec t couplings. Finally 
we put the constant factor of TV — 2 into the equation to adjust the energy back 
to zero overall for valid solutions which contain N — 2 connections. 

3.2.3 Construction of E over i ap (q) 

It is much easier to prevent amino acids from overlapping in the diamond map- 
ping than in the turn mapping. The only way that amino acids could overlap 
in the diamond mapping is for amino acids which have an even number of 
bonds between them to flip bits corresponding to the same lattice location. For 
instance, in Fig. 14 its clear that the fourth amino acid could overlap with sec- 
ond amino acid since the orange lattice points are possibilities for both amino 
acids. Assuming that the diamond lattice positions are encoded so that the 
inner diamond bits come first in the bit-string for each amino acid and that bits 
are enumerated in some consistent fashion (e.g. starting at the top and going 
clockwise around the diamond), we can write the following, 

N-l N n k 

Kverlap (<?) = Kverlap ^ EE K 1 + k ~ ^ m ° d ^ ^ ' ^ 

k=2 h>k i=l 

This expression would perfectly sum over all the possible overlaps as the first two 
sums iterate through all possible overlapping pairs and the third sum iterates 
through all of the diamond points up to the last point they both share, n^. 

3.2.4 Construction of E pair (q) 

To form the pairwise interaction term we simply couple each bit to the possible 
adjacent lattice locations which could be occupied by other amino acids. The 
strength of the coupling will depend on the interaction matrix element between 
the two amino acids coupled by the term. Additionally, we note that amino 
acids are only able to be adjacent if there are an even number of amino acids 
(2 or greater) in between the two. Thus, the formula is as follows: 

N-l N 

E pair (<?) = E E E J ^ P - h ) mod 2 ] ihj (46) 

k=2 h=k+2 <ij> 

where the sum over < ij > is understood as a sum over bits corresponding to 
adjacent lattice sites. There is no straightforward way to write the function 
< ij > in analytical terms. Nevertheless, for large problems it is trivial to write 
a program which iterates through bits in the second amino acid with a for-loop 
and evaluated the sum on those bits if the first amino acid bit and the second 
amino acid bit have a grid distance of 1. 
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4 Pseudo-boolean Function to W-SAT 



In order to take advantage of state-of-the-art satisfiability (SAT) solvers to 
optimize our pseudo-boolean function, it is necessary to map the problem to 
Weighted Maximum Satisfiability (W-SAT). The most general form of the generic 
SAT problem is known as K-SAT. In K-SAT the problem is to find a vector of 
boolean valued variables which satisfies a list of clauses, each containing up to 
K variables, which constrain the solution. When K-SAT has a solution it is 
known as "satisfiable" and for K < 2 the problem is tractable in polynomial 
time. However, for K > 2 the problem is known to be NP-complete; in fact, 
3-SAT was the first problem proved to be NP-Complete (Cook 1971). 



Maximum Satisfiability (MAX-SAT) is a more difficult version of the canonical 
SAT problem which is relevant when K-SAT is either "unsatisfiable" or at least 
not known to be satisfiable. In MAX-SAT the goal is not necessarily to find the 
solution string which satisfies all clauses (such a solution string may not even 
exist); rather, the goal is to find the solution string which satisfies the maximum 
number of clauses. 

An extension of MAX-SAT known as Weighted Maximum Satisfiability (aka 
W-SAT) is what will be most relevant to us. In W-SAT each clause is given 
a positive integer valued "weight" which is added to a sum only if the clause 
evaluates to false. Accordingly, in W-SAT the goal is to minimize this sum 
rather than the total number of false clauses as in canonical MAX-SAT (Xing 
& Zhang 2005, Boros & Hammer 2002). We can more succinctly state the 
problem as follows: given m number of clauses (y) each with a weight of w, 
minimize 



The same approximation schemes and exact solver algorithms which work 
for MAX-SAT also work for W-SAT (Boughaci & Drias 2004, Pankratov & 
Borodin 2010). In order to use these solvers one must first translate their 
pseudo-boolean function into a W-SAT problem articulated in what is known as 
Weighted Conjunctive Normal Form (WCNF). In WCNF, the W-SAT problem 
is stated as a list of weights followed by a clause with each clause stated as 
an OR statement between integers representing the index of the corresponding 
boolean variable in the solution vector. In WCNF, a negative integer denotes a 
negation. For instance the WCNF clause "4000 9 —1 82" means igV^iVz^ 
with penalty of 4000 if clause evaluates to false. Fig. 15 shows this clause as 
a logic circuit. 



4.1 MAX-SAT and W-SAT 




(47) 
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(1 - Xl + XXXS2 + XxXg - XxXCjX^) 



X&2 



Figure 15: A logical circuit representation of the CNF clause: "9 —1 82" 
4.2 Constructing WCNF clauses 

To prepare the WCNF input file from a pseudo-boolean function one will need to 
write a short script which transforms each term in the pseudo-boolean function 
into a WCNF clause. There is more than one way to accomplish this transfor- 
mation and we will only discuss one method here. For a more complete review 
of this topic, see (Een & Sorensson 2006). 

It will be very useful to think of CNF clauses as logical circuits which involve 
only OR gates and not gates as in Fig. 15. Weights in WCNF notation always 
represent a positive value. Because pseudo-boolean functions are treated as cost 
functions to minimize and the goal of W-SAT is to minimize the sum of weights 
on false clauses, terms in the pseudo-boolean function with a positive weight 
are very easy to translate in WCNF notation. To achieve this, one needs only 
to pass all variables in the clause through a not gate and then a series of OR 
gates (effectively making a NAND gate which takes all variables as input). This 
circuit is illustrated in Fig. 16 for the case of a 5 variables clause. 




(1 - xix 2 xzXiX 5 ) 



Figure 16: A logical circuit which shows that any pseudo-boolean term with 
positive weight is equivalent (up to a constant) to a CNF clause with each 
variable negated. The term produced here is negative because the weight is 
only added when the clause evaluates to false. 

Representing a negative weighted pseudo-boolean term in CNF is less trivial 
but follows a simple pattern. To make the CNF clause positive (corresponding 
to negative boolean term) one needs to construct the same circuit as in the case 
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when the boolean term is positive but remove one of the NOT gates. An example 
comprising three variables is shown in Fig. 17. However, this circuit alone does 




(1 - X X X 2 + X 1 X 2 X 3 ) 



Figure 17: A logical circuit on three variables which gives a positive valued 
3-local CNF term. 

not accomplish our goal as it produces a 2-local term with negative weight in 
addition to the 3-local term with positive weight. Consequentially, after using 
the circuit in Fig. 17 to get rid of the 3-local term ll x\x 2 x^ we must subtract 
the term "£1X2" multiplied by its weight from the pseudo-boolean expression 
we are converting into CNF. At first glance, it is not obvious that this procedure 
will get us anywhere - we turned a term into CNF only to introduce a new term 
into the pseudo-boolean which we must convert back into CNF. However, the 
auxiliary terms produced by this circuit are of one degree less than number of 
variables in the term; thus, we can iterate this procedure until only the constant 
term remains. The next CNF clause (this time 2-local) is shown in Fig. 18. 



Figure 18: A logical circuit on three variables which gives a positive valued 
2-local CNF term. 

4.3 Solving SAT problems 

While MAX-SAT is known to be NP-hard, there exist heuristic algorithms which 
are guaranteed to satisfy a fixed fraction of the clauses of the optimal solution 
in polynomial time. In general, oblivious local search will achieve at least an 
approximation ratio of Tabu search achieves a ratio of at least |±| and 

non-oblivious local search achieves an approximation ratio of 2 2 ^T 1 where k 
is the "K" in K-SAT. For the special case of MAX-2-SAT the best possible 
algorithm is theoretically capable of satisfying at least §5 + e ~ 0.955 + e [9] in 
polynomial time (Pankratov & Borodin 2010, Choi, Standley & Darwichc 2009). 
Additionally, there are a great deal of exact MAX-SAT solvers which run in 
super-polynomial time but in many cases can find the solution to MAX-SAT in 
a very short amount of time, even for problems containing hundreds of variables 
and clauses (Marques-Silva & Sakallah 2007, Larrosa, Heras & De Givry 2006). 
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5 W-SAT to Integer-Linear Programming 



Integer-Linear Programming (ILP) is a subset of linear programming problems 
in which some variables are restricted to integer domains. In general, ILP is an 
NP-Hard problem but the importance of ILP problems (particular for logistics 
scheduling) has produced many extremely good exponential-time exact solvers 
and polynomial-time heuristic solvers (Xing & Zhang 2005). Pseudo-boolean 
optimization is an even more specific case of ILP sometimes known as 0-1 ILP 
where the integer variables are boolean (Boros & Hammer 2002) . The mapping 
between W-SAT and ILP is very straightforward. 

5.1 Mapping to ILP 

In ILP, the goal is to minimize an objective function of integer-valued variables 
subject to a list of inequality constraints which must be satisfied. The inequality 
constraints come directly from the clauses in W-SAT. As described in Sec. 4.1, 
the logical clause from the WCNF clause "4000 9 -1 82" (which again, means 
xgV^xi VX82 with penalty of 4000 if clause evaluates to False) can be represented 
as x g + (1 — xi) + x S 2 > 1 s.t. x n e {0, 1}. In ILP, all constraints must be 
satisfied but in W-SAT clauses are sometimes not satisfied; to accommodate 
this we introduce an auxiliary binary variable, 2/1 into the equation and get 
2/i + xg + (1 — Xi) + x S 2 > 1- Thus, if the original equation is False, 2/1 will 
have a value of True which satisfies the inequality. We can take advantage of 
this auxiliary variable to construct the optimization function, W. Since the 
clause in our example has a weight of 4000 we can write W = 4000j/i s.t. 
j/i + x 9 + (1 — xi) + x%2 > 1- Thus, the mapping between ILP and W-SAT 
is extremely trivial: all WCNF clauses are rewritten as linear equalities which 
are > 1 — 2/i by adding together the variables (or their negations) where i is 
the index of the clause and the objective function is written as W — J2iLi w iVi 
where N is the number of clauses and Wi is the weight of that clause (Xing & 
Zhang 2005). 

5.2 Solving ILP problems 

Commercial logistic scheduling software such as IBM ILOG CPLEX Optimiza- 
tion Studio (aka CPLEX) is designed to solve in integer programming, linear 
programming, and mixed integer- linear programming problems on a very large 
scale (IBM 2009). Constraint satisfaction problems which are sometimes very 
difficult to solve using conventional SAT techniques can be easier to solve us- 
ing ILP techniques and vice versa. In particular, SAT solvers and specialized 
pseudo-boolean optimizers seem to outperform ILP solvers when a problem 
is over-constrained (Aloul, Ramani, Markov & Sakallah 2002). On the other 
hand, for problems which are under-constrained and have a large number of 
variables ILP solvers are the natural choice. In some cases 0-1 ILP optimizers 
such as Pueblo will outperform both SAT solvers and commercial ILP solvers 
(Sakallah 2006, Sheini & Sakallah 2005, Manolios & Papavasileiou 2011). 
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6 Locality Reductions 



The practical ability to either exactly or approximately solve random instances 
of constraint satisfaction optimization such as pseudo-boolean optimization or 
MAX-SAT seems to depend very sensitively on the variable to clause ratio and 
degree of constraint expressions (Pankratov & Borodin 2010, Kahl & Strandmark 
2011, Xing & Zhang 2005). In fact, the degree of constraints determines the com- 
plexity class of certain constraint satisfaction problems; e.g. 2-SAT is proven 
to be in P whereas 3-SAT is in NP-Complcte (Cook 1971). Clearly for in- 
stances such as this there can be no efficient method which reduces the degree of 
constraints. Fortunately, reducing the degree of constraints in general pseudo- 
boolean optimization (i.e. reducing the polynomial order of pseudo-boolean 
terms) can be done efficiently. 

Constraint degree reduction is particularly important if we wish to solve 
our problem using existing architectures for adiabatic quantum computation 
because available devices tend to be very limited in their ability to realize arbi- 
trary variable couplings (especially high ordered couplings). For instance, the 
D-Wave One device used for pseudo-boolean optimization in (Perdomo-Ortiz 
et al. 2012) is only able to implement 2-local qubit couplings and has limited 
coupler resolution. To encode functions of higher locality in such setups, we 
must introduce ancilla bits which replace 2-local terms to reduce locality. Be- 
cause these ancilla become free parameters of the system, it is also necessary to 
introduce penalty functions to account for the possibility that their value may 
be incorrect. All of this is accomplished with the function E A (qi,qj,q n ',d n ) in 
Eq. 48 which introduces the ancillary bit q n in order to collapse the 2-local 
term q^j with energy penalty of S n if q n ^ qiqj. For a further discussion, see 
(Perdomo et al. 2008). 

E A (qi, q 3 ,q n ; S n ) = S n (3q n + q t qj - 2q { q n - 2q j q n ) (48) 

If one desires an entirely 2-local energy function then many E A (<?,, qj,q n ; <5„)'s 
may be necessary to collapse all high-local terms. For instance, consider the 
complete energy function for the HP model protein HPPHP when coded in 
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the turn ancilla mapping: 
E = —4(7296^1 + 4gig 3 g 6 Ai + 3geAi + 28giA 2 + 25gig2A 2 + 108g2A 2 — 56qiqs\2 (49) 

— 50gig 2 g 3 A2 + 26g 2 g 3 A 2 + 28q 3 \ 2 + 24gig 4 A 2 — 16gig 2 g4A 2 — 56g 2 g4A 2 — 48gig 3 g 4 A 2 
+ 32gig 2 g 3 g 4 A 2 — 18g 2 g 3 g 4 A 2 + 25g 3 g 4 A 2 + 108g 4 A 2 — 56q 1 q 5 X 2 — 48qiq 2 q 5 \2 

+ 25q 2 q 5 \ 2 + 48gig 3 g 5 A 2 - 50q 2 q 3 q 5 \ 2 - 56g 3 gsA 2 - 48gig 4 g 5 A 2 + S2qiq 2 q 4 qsX2 

— 18q 2 q4qr,\ 2 + 36g 2 g 3 g4g5 A 2 — 50g 3 g4gsA 2 + 25g4gsA 2 + 28goA 2 — 32gig 7 A 2 

— 96g 2 g7A 2 + 64gig 3 g7A 2 — 32g 3 g7A 2 + 64g 2 g4g7A 2 — 96g4g7A 2 + 64gigsg7A 2 

+ 64g 3 g 5 g 7 A2 — 32g 5 g 7 A2 — 32g 7 A2 — 16g±g 8 A 2 — 48g2gsA2 + 32gig 3 g 8 A 2 — 16g 3 ggA2 
+ 32g 2 g4g8A 2 - 48g 4 ggA 2 + 32gig 5 g 8 A 2 + 32g 3 g 5 g8A 2 - 16g 5 g8A 2 + 64g 7 ggA2 

— 32g 8 A 2 — 8gig;)A 2 — 24g 2 g9A2 + 16gig 3 ggA2 — 8g 3 g9A 2 + 16g2g 4 ggA2 — 24g 4 ggA2 
+ 16gig 5 g9A 2 + 16g 3 g 5 ggA 2 — 8g 5 ggA 2 + 32g 7 ggA 2 + 16g 8 g9A 2 — 20ggA 2 — 4gigi A 2 

— 12g 2 gi A 2 + 8gig 3 gi ( )A 2 - 4g 3 gioA 2 + 8g 2 94gioA 2 - 12g 4 gioA 2 + 8gig 5 gi A 2 
+ 8g 3 g5gioA 2 — 4g 5 gi A 2 + 16g7gioA 2 + 8g 8 gioA 2 + 4g 9 gi A 2 — llgioA 2 + 36A 2 . 

In order to reduce this function to 2-local we will need to collapse some of 
the 2-local terms inside of the 3- local terms to a single bit. We enumerate all of 
the 3-local terms and their corresponding 2-local terms which we could use to 
reduce each 3-local term in Eq. 50. 
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Eq. 50 shows that there are 30, 3-local terms in Eq. 49 and three different 
ways to collapse each of those 3-local terms. In general, the problem of choosing 
the most efficient 2-local terms to collapse this function is NP-Complete. This 
becomes evident if we represent our problem as an element cover on a bipartite 
graph. Suppose we relabel each 3-local term on the left as "set" 1-30, denoted 
as S\S<2---Szq. We can then make the following bipartite graph which connects 
the 3-local terms to the 2-local terms which collapse them. 

Fig. 21 shows that we can now restate the problem in the following way: 
"choose the fewest number of 2-local terms (on the left) which covers all 3-local 
terms (on the right) with at least one edge." In general, this problem is isomor- 
phic to the canonical "hitting set" problem which is equivalent to set cover, one 
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Figure 19: A bipartite graph connecting the 3-local terms (S n ) in Eq. 49 to the 
2-local terms {qiqj) which collapse them. 

of Karp's 21 NP-Complete Problems (Chandrasekaran, Karp, Moreno- Centeno 
& Vempala 2011, Chvatal 1979, Laue 2008). However, we have specifically kept 
this issue in mind when creating the turn-ancilla representation in such a way 
as to guarantee that it is easy to find a relatively efficient solution to this prob- 
lem. Accordingly, our experience has been that a greedy local-search algorithm 
performs very well. 

The explanation for this is simple: each 3 or 4-local term will contain no 
more than 1 ancillary bit; thus, to cover all 3 and 4-local terms we can focus 
entirely on the physical bits (in this case, bits 1-5). In alternative mappings not 
presented here we have frequently encountered extremely difficult instances of 
the hitting set problem during the reduction process. In these situations one 
should see (Shi & Cai 2010) for a very efficient algorithm which can exactly 
solve hitting cover in O (1.23801"). 
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7 Quantum Realization 



A primary goal of this review is to elucidate an efficient process for encoding 
chemical physics problems into a form suitable for quantum computation. In 
addition to providing the alternatives for the solution of the lattice heteropoly- 
mer problem in quantum devices, we seek to provide a general explanation of 
considerations for constructing energy functions for these devices. These have 
many possible applications for solving problems related to statistical mechanics 
on the device. In this section, we will complete our review by demonstrating 
the final steps required to embed a small instance of a particular lattice protein 
problem into a QUBO Hamiltonian. 

The Hamiltonians and the number of resources presented here correspond 
to the minimum amount of resources needed assuming the device can handle 
many-body interactions as is the case for NMR quantum computers or trapped 
ions. The hierarchical experimental proposals presented here work for lattice 
folding under no external constraints, i.e., amino acid chains in "free space" 1 . 
As a final step we will reduce these Hamiltonians to a 2-local form specifically 
design for the Dwave One used in (Pcrdomo-Ortiz et al. 2012, Neven, Rose & 
Macready 2008, Denchev, Ding, Vishwanathan & Neven 2012, Johnson, Amin, 
Gildcrt, Lanting, Hamzc, Dickson, Harris, Berkley, Johansson, Bunyk, Chappie, 
Endcrud, Hilton, Karimi, Ladizinsky, Ladizinsky, Oh, Pcrminov, Rich, Thorn, 
Tolkacheva, Truncik, Uchaikin, Wang, Wilson & Rose 2011). The final Hamilto- 
nian we present is more efficient than that used in (Perdomo-Ortiz et al. 2012) as 
we have since realized several tricks to make the energy function more compact. 

7.1 Previous experimental implementation 

Throughout this review we have referred to an experimental implementation of 
quantum annealing to solve lattice heteropolymer problems in (Perdomo-Ortiz 
et al. 2012). The quantum hardware employed consists of 16 units of a re- 
cently characterized eight qubit unit cell (Johnson et al. 2011, Harris, Johnson, 
Lanting, Berkley, Johansson, Bunyk, Tolkacheva, Ladizinsky, Ladizinsky, Oh, 
Cioata, Perminov, Spear, Enderud, Rich, Uchaikin, Thorn, Chappie, Wang, Wil- 
son, Amin, Dickson, Karimi, Macready, Truncik & Rose 2010). Post-fabrication 
characterization determined that only 115 qubits out of the 128 qubit array can 
be reliably used for computation. The array of coupled superconducting flux 
qubits is, effectively, an artificial Ising spin system with programmable spin- 
spin couplings and transverse magnetic fields. It is designed to solve instances 
of the following (NP-hard) classical optimization problem: given a set of lo- 
cal longitudinal fields (hi) and an interaction matrix (Jij), find the assignment 
s = S!S2S 3 ...s N , that minimizes the objective function E(s), where, 

E(s) = ^ M» + ^ JijSiSj (51) 

l<i<N l<i<j<N 

1 External interactions could also be included as presented and verified experimentally in 
(Perdomo-Ortiz ct al. 2012). 
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and Si € —1,1. Thus, the solution to this problem, s, can be encoded into the 
ground-state wavefunction of the quantum Hamiltonian, 

l<i<W l<i<J<JV 

Quantum annealing exploits the adiabatic theorem of quantum mechanics, 
which states that a quantum system initialized in the ground state of a time- 
dependent Hamiltonian remains in the instantaneous ground state, as long as it 
is driven sufficiently slowly. Since the ground state of H p encodes the solution to 
the optimization problem, the idea behind quantum annealing is to adiabatically 
prepare this ground state by initializing the quantum system in some easy-to- 
prepare ground state, Hb- In this case, Hb corresponds to a superposition of all 
states of the computational basis. The system is driven slowly to the problem 
Hamiltonian, %{t = 1) "Hp. Deviations from the ground-state are expected 
due to deviations from adiabaticity, as well as thermal noise and imperfections 
in the implementation of the Hamiltonian. 

Using the encoding methods discussed here, the authors were able to encode 
and to solve the global minima solution for small tetrapeptide and hexapeptide 
chains under several experimental schemes involving 5 and 8 qubits for four- 
amino-acid sequence (Hydrophobic-Polar model) and 5, 27, 28, and 81 qubits ex- 
periments for the six-amino-acid sequence under the Miyazawa-Jernigan model 
for general pairwise interactions. 



7.2 Six unit Miyazawa-Jernigan protein 

The example we will present here is a different encoding of the largest problem 
performed in (Perdomo-Ortiz et al. 2012): the Miyazawa-Jernigan (MJ) protein, 
Proline-Serine- Valine-Lysine-Methionine- Alanine (PSVKMA) on a 2D lattice. 
We will use the pair-wise nearest-neighbor MJ interaction energies presented in 
Table 3 of (Miyazawa & Jernigan 1996) and shown in Fig. 20. 
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Figure 20: Interaction matrix for our protein in the MJ model. 



We will use the turn ancilla construction for our energy function and constrain 
the first three virtual bits to 010, as before. Recall that the turn ancilla con- 
struction requires 2N — 5 physical information bits; thus, our 6-unit MJ protein 
will be encoded into 7 bits. 
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7.2.1 E back {q) for 6-unit SAW on 2D lattice 

Using Eq. 13, we find that our 6-unit protein has the backwards energy function, 

Eback(q) = hack(qiq2 - 2(71(7392 + 2q 3 q 2 - 2q 3 q 4 q 2 - 2q 3 q 5 q 2 (53) 
+ 4^394(7592 - 2949592 + 9592 + 9394 - 2(7 3 <74 95 + 29495 - 2949596 
+ 9596 + 9497 - 2949597 - 2949697 + 494959697 - 2959697 + 9697)- 

Soon, we will discuss how to choose the appropriate value for \b ac k but for now 
we simply note that Xback and X veria P penalize the same illegal folds; thus we 
realize that Xback — X over i a p. 

7.2.2 E oveHap (q) for 6-unit SAW on 2D lattice 

Using Eq. 30, we calculate the overlap energy function as, 

Eoveriap(q) = X over ia P (96<?2 qi — 96(j2g39i — 64^3^1 — 64g 2 g4gi + 64529354(71 — 96<73949i + 96g4(?i 

- 96<72(?59i + 64g 2 (74959i — 96<74(75<7i — 64g 5 gi — 48Q29691 + 32g 2 g 3 g 6 gi — 48g 3 g 6 gi + 32g 3 g 4 g6gi 

- 48g4?69i + 32g 2 95969i + 32g 4 gsg69i — 48gs969i + 72g 6 gi — 48g2g7gi — 48g3g7gi + 32g 2 g4g7gi 

- 48g 4 g7gi + 96g 3 gsg7gi - 48g 5 g7gi + 32g 2 g 6 g7gi + 32g 4 geg7gi - 48g 6 g7gi - 8g 7 gi - 8q 3 qw 

+ 64g 3 g 8 gi + 64g 5 g 8 gi - 32g 8 gi + 32g 3 g 9 gi + 32g 5 g 9 gi - 16g 9 gi + 16g 3 giogi + 16g 5 giogi - 8gi gi 

+ 8g3gngi + 8g 5 giigi — 4gngi + 64g 3 gi2gi + 64g 5 gi 2 gi + 64g 7 gi2gi — 96gi 2 gi + 32g 3 gi 3 gi 

+ 32g 5 gi 3 gi + 32g 7 gi 3 gi - 48gi 3 gi + 16g 3 gi 4 gi + 16g 5 gi4gi + 16g 7 gi 4 gi - 24g i4 gi + 8q 3 q 15 qi 

+ 8g59i5gi + 8g 7 gi5gi - 12gi 5 gi + 64gi + 144g 2 + 96g 2 g 3 + 64g 3 - 64g 2 g 4 - 64g 2 g3g4 + 96g 3 g 4 + 144g 4 

+ 96g 2 g5 - 96g 2 g3gs - 64g 3 g 5 - 64g 2 g4gs + 64g 2 g3g4gr, - 96g 3 g4gs + 96g 4 gs + 64g 5 - 8g 2 g 6 - 48g 2 g3g6 

+ 72g 3 g 6 - 48g 2 g4ge - 48g 3 g4g6 - 8g 4 ge - 48g 2 g 5 g6 + 32g 2 g3g5g6 - 48g 3 g 5 g6 + 32g 3 g 4 gsg6 - 48g 4 gs96 

+ 72g 5 g 6 + 36g 6 + 72g 2 g7 - 48g 2 g3g7 - 8g 3 g7 - 48g 2 g4g7 + 32g 2 g 3 g 4 g7 - 48g 3 g4g7 + 72g 4 g7 - 48g 2 gsg7 

- 48g 3 g5g7 + 32g 2 g4gsg7 - 48g 4 gsg7 - 89597 - 48g 2 g 6 g7 + 32g 2 g396g7 - 48g 3 g6g7 + 32g 3 g 4 g6g7 
+ 32g 2 g5g6g7 + 32g 4 g5geg7 - 48g 5 g6g7 + 72g 6 g7 + 36g 7 - 96g 2 g 8 - 32g 3 g 8 + 64g 2 g 4 g8 - 96g 4 g 8 

- 32g 5 g 8 — 32g 8 — 48g2g9 — 16g 3 gg + 32g 2 g 4 g9 — 48g4g 9 + 32g 3 g 5 g9 — 16g 5 g9 + 64g 8 g 9 — 32g 9 — 24g 2 gio 
+ 16g 2 g4gio - 24g 4 gi + 16g 3 g 5 gi - 8g 5 gi + 32g 8 g 10 + 16g 9 gi - 20gi - 12g 2 gu - 4g 3 gn + 8g 2 g4gn 

- 12q 4 qn + 8g 3 g 5 gn - 4g 5 gn + 16g 8 gn + 8g 9 gn + 4gi gn - llgn - 96g 2 gi 2 - 96g 3 gi 2 + 64g 2 g 4 gi 2 

- 96g 4 gi2 + 64g 3 g 5 gi 2 - 96g 5 gi2 + 64g 2 g6gi2 + 64g 4 g6gi2 - 96g 6 gi2 + 64g 3 g 7 gi2 + 64g 5 g 7 gi2 

- 96g 7 gi2 + 64g i2 - 48g 2 gi 3 - 48g 3 gi 3 + 32g 2 g4gi 3 - 48g 4 gi3 + 32g 3 g 5 gi 3 - 48g 5 gi3 + 32g 2 g6gi3 
+ 32g 4 g 6 gi3 - 48g 6 gi 3 + 32g 3 g 7 gi 3 + 32g 5 g 7 gi 3 - 48g 7 gi3 + 64gi 2 gi 3 + 16gi 3 - 24g 2 gi4 - 24g 3 gi 4 
+ 16g2g4gi4 — 24g 4 gi4 + 16g3gsgi4 — 24g 5 gi 4 + 16g2g6gi4 + 16g4g6gi4 — 24g 6 gi 4 + 16g3g7gi4 

+ 16g 5 g7gi4 - 24g 7 gi4 + 32gi 2 gi 4 + 16gi 3 gi4 + 4g i4 - 12g 2 gi 5 - 12g 3 gi 5 + 8g 2 g4gi5 - 12g 4 gis 
+ 8g 3 g5gi5 - 12g 5 gi5 + 8g 2 g6gis + 8g 4 g6gis - 12g6gis + 8g 3 g7gis + 8g 5 g7gis - 12g 7 gi5 
+ 16gi 2 gi5 + 8gi 3 gi5 + 4gi 4 gi 5 + gis - 48g 4 geg7 + 64g 3 g5g 8 ). (54) 

We notice that as discussed in Sec. 6, all the 3-local terms here contain at least 
two physical information qubits (i.e. qi through qi). 
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7.2.3 E pair (q) for MJ-model PSVKMA 

Using the J matrix as denned in Eq. 32 we calculate the pair-wise energy func- 
tion as, 

Ep a ir(q) = — 4q 2 gi6 + 4gi<739i6 + 3gi6 — 8171517 — 1692917 + 89193917 — 853917 

+ 89294917 — I694917 + 89195917 + 89395917 — 895917 + 89296917 + 89496917 — 1695917 

+ 89197917 + 89397917 + 89597917 - 897917 + 30gi7 - 129191s - 12g 2 gi 8 + 12919391s 

- 1293918 + 129294918 - I294918 + 12gi9 5 gi 8 + 12g3g 5 gi 8 - 12g 5 gi 8 + 21gi 8 - 16g 2 gi9 

- I693919 + I69294919 - I694919 + I69395919 - I695919 + I69296919 + I69496919 

- 16g 6 9i9 + I69397919 + I69597919 - I697919 + 28919. (55) 

7.2.4 Setting A penalty values 

Finally, we will discuss how one chooses the correct penalty values for the en- 
ergy function. This is a crucial step if one wishes to implement the algorithm 
experimentally as all currently available architectures for adiabatic quantum an- 
nealing have limited coupler resolution. That is, quantum annealing machines 
cannot realize arbitrary constant values for the QUBO expression. Thus, it is 
very important that one chooses the lowest possible penalty values which still 
impose the correct constraints. In our problem we choose the value of X OV eriap 
by asking ourselves: what is the greatest possible amount that any overlap could 
lower the system energy? In general, a very conservative upper bound can be 
obtained by simply summing together every J matrix element (which would 
mean that a single overlap allowed every single possible interaction to occur); in 
our problem this upper-bound would be -10. Thus, we can set \ OV eria P = +10. 



7.2.5 Reduction to 2-local 

Using a standard greedy search algorithm we find that an efficient way to col- 
lapse this energy function to 2-local is to make ancilla with the qubit pairs, 

9294 -> 920 

9i93 -> 921 

9395 922 
91 95 -> 923 

9296 -> 924 (56) 

9496 -> 925 
9397 -> 926 
9597 -> 927 
91 97 928 • 

There is one issue left to discuss - the value of 8 n in Eq. 48. The purpose of 
S n is to constrain the reductions in Eq. 56 so that the value of the ancillary bit 
actually corresponds to the product of the two bits it is supposed to represent. 
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Table 1: Truth table for the function E A (qi, qj,q n ; S n ) from Eq. 48. 

q n qi q 3 E A (q u q 3 ■, q n ; 5 n ) 



1 

10 
111 

1 3(5 

10 1 5 
110 5 

011 8 



To understand how Eq. 48 accomplishes this see Table 7.2.5. In order for Eq. 48 
to work we must choose S n which is large enough so that a violation of the 
reduction we desire will always raise the system energy. Thus, we must ensure 
that S n is large enough so that configurations which do not conform to the 
reduction are penalized by an amount higher than the largest penalty they 
could avoid and larger in magnitude than the largest energy reduction they 
could achieve with the illegal move. Of course, finding the exact minimum value 
of E{q) is as difficult as minimizing E{q) (our goal). Instead, we can simply 
make an upper-bound for the penalty by setting it equal to one plus either the 
sum of the absolute value of all negative psuedo-boolean coefficients or the sum 
of all positive psuedo-boolean coefficients corresponding to the variables being 
collapsed in E{q) (whichever sum is larger). 

7.2.6 QUBO Matrix and Solutions 

After reduction of the energy function to 2-local, we arrive at the final pseudo- 
boolean energy function. Instead of writing out the entire pseudo-boolean ex- 
pression we will instead provide a matrix containing all of the coefficients of 
1-local terms on the diagonal and 2-local terms in the upper triangular portion 
of this matrix. This representation is known as the QUBO matrix and con- 
tains all of the couplings needed for experimental implementation and is shown 
in Eq. 57. Note that the full pseudo-boolean expression contains one constant 
term that we drop in the matrix representation. This constant has a value of 
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C = 180 for this particular problem. 
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(57) 

Taking the matrix in Eq. 57 as Q, we can write the total energy of a given 
solution (denoted by q) as, 

E(q) = qQq. (58) 

The problem is now ready for its implementation on a quantum device. For our 
particular problem instance the solution string is given by the bit string, 

0,0,0,1,0,1,1,1,1,0,0,1,1,1,0,0,1,0,1,0,0,0,0,0,1,0,0,0. (59) 

The energy given by Eq. 58 is —186. In the original expression this corresponds 
to an energy of C — 186 = 180 — 186 = —6. Let's confirm that this is accurate to 
the MJ model. Looking only at the physical information bits and prepending the 
first three constant bits (010) we see that the bit string prescribes the following 
fold: 

q = J^^J^^^ (60) 

right down down left up 

which corresponds to the fold, 




Figure 21: The solution to our example problem for MJ protein PSVKMA. 
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8 Conclusion 



As both traditional and quantum computer science continue to advance as fields, 
domain scientists from all disciplines need to develop new ways of representing 
problems in order to leverage state-of-the-art computational tools. In this re- 
view, we discussed strategies and techniques for solving lattice heteropolymer 
problems with some of these tools. While the lattice heteropolymer model is 
widely applicable to many problems, the general principles used to optimally 
encode and constrain this particular application are fairly universal for discrete 
optimization problems in the physical sciences. 

We focused on three mappings: "turn ancilla" , "turn circuit" and "dia- 
mond" . The turn ancilla mapping is the best mapping in terms of the scaling of 
the number of resources for large instances, thus making it ideal for benchmark 
studies of lattice folding using (heuristic) solvers for pseudo-boolean minimiza- 
tion. Additionally, this method shows how one can use ancilla variables to 
construct a fitness function with relatively few constraints per clause (i.e. low- 
locality). With ancilla variables even an extremely simple encoding, such as 
the turn encoding, can be used to construct a complicated energy function. 
While some of the particular tricks employed to optimize the efficiency of this 
mapping, such as introducing the backwards penalty, are specific to lattice het- 
eropolymcrs, the general logic behind these tricks is much more universal. 

The turn circuit mapping is the most compact of all three mappings. The 
extremely efficient use of variables (qubits) makes it ideal for benchmark exper- 
iments on quantum devices which can handle many body couplings. Moreover, 
the turn circuit method demonstrates how one can construct an elaborate energy 
function by utilizing logic circuits to put together a high-local fitness function 
of arbitrary complexity without ancilla variables. While different problems may 
involve different circuits, the underlying strategy is very broadly applicable. 

The diamond encoding illustrates a strategy for producing an extremely 
under-constrained optimization problem. Furthermore, this method demon- 
strates that even fairly complex energy functions can be represented as natively 
2-local functions if one is willing to sacrifice efficiency. Many quantum devices 
can only couple bits pairwise; thus, this is a very important quality of the dia- 
mond encoding. Finally, if one uses another, more efficient encoding, we explain 
how reductions can be used to replace high-local terms with 2-local terms in an 
optimally efficient fashion but at the cost of needing very high coupler resolu- 
tion. The relatively few constraints in the diamond encoding make it a natural 
choice for exact or heuristic ILP and W-SAT solvers. 

These three strategies elucidate many of the concepts that we find important 
when producing problems suitable for the D-Wave device utilized in (Perdomo- 
Ortiz ct al. 2012). Accordingly, as quantum information science continues to 
develop, we hope that the methods discussed in this review will be useful to 
scientists wishing to leverage similar technology for the solution of discrete op- 
timization problems. 
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